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Abstract. We report here the first full sky component separation and CMB power spectrum estimation using a 
Wiener filtering technique on simulated data from the upcoming MAP experiment, set to launch in early 2001. 
The simulations included contributions from the three dominant astrophysical components expected in the five 
MAP spectral bands, namely CMB radiation, Galactic dust, and synchrotron emission. We assumed a simple 
homogeneous and isotropic white noise model and performed our analysis up to a spherical harmonic multipole 
^max = 512 on the fraction of the sky defined by |6| > 20°. We find that the reconstruction errors are reasonably 
well fitted by a Gaussian with an rms of 24 fiK, but with significant deviations in the tails. Our results further 
suppo rt the predictions on the resulting CMB power spectrum of a previous estimate by Bouchet & Gispert 



199£), which entailed a number of assumptions this work removes. 
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1. Introduction 

The upcoming CMB satellite experiments MAP and the 
Planck Surveyor offer an unprecedented opportunity to 
measure CMB temperature fluctuations on the whole 
sky. A major hurdle in extracting the primary CMB sig- 
nal from data, apart from noise, is the removal of the 
Galactic and extragalactic foregrounds. However, as the 
foregrounds differ from CMB in both frequency depen- 
dence and spatial distribution, one can reduce their resid- 
ual level in a multi- frequency CMB experiment. A multi- 
frequency, multi-resolution, Wiener filtering method to 
optimally extract the CMB component was developed 



(Bouchet et al. 1996; Tegmark & Efstathiou 1996; Bouchet 



fe Gispert 199£) and lead to detailed predictions for the 
accuracy achievable with several experiments, and in par- 
ticular, for MAP QBouchct fc Gispert 1999| ). These pa- 
pers were purely semi-analytical and entailed a number of 
simplifying approximations and assumptions. A practical 
implementation was performed by Bouchet et al. (1996) 
on the simulated data produced by Bouchet et al. (1995| ). 



These numerical studies confirmed that the residual con- 
tamination after cleaning the map is much smaller than 
the CMB primary signal and therefore the foregrounds 
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may not be a major obstacle in the extraction of CMB 
temperature angular power spectrum if the real sky be- 
haves similarly to the assumed model. 

These simulation studies were based on filtering the 
data obtained on small patches of the sky of typically 10° 
by 10° angular size. They did not involve the assumption 
of Gaussian foregrounds since actual templates were used 
(e.g. portions of the Haslam and IRAS maps at respec- 
tively 408 GHz and 100 [im). The results, however, could 
not be safely extrapolated to the full sky as a result of 
three facts: 1) the Wiener filter used was always optimal 
for the particular (small) region being analyzed instead of 
being an average filter for all the sky; 2) nothing could be 
said about modes larger than the map size; 3) there was 
no really satisfactory way of co-adding individual map re- 
sults even if all maps covering the sky had been used (one 
problem being the use of Fourier transforms with peri- 
odic boundary conditions). Note that a similar approach 
based on a multi-frequency maximum entropy method in 
Fourier space on the same simulations leads to compara- 
ble optimistic conclusions ( Hobson et al. (1998] ), applied 
to simulated MAP data by lones et al. (1999^ but with 
similar caveats. 

In this work, we extend the Wiener filtering method 
to full-sky maps with a symmetric Galaxy cut (\b\ > 20°) 
and discuss our solutions to the various technical difficul- 
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ties encountered. We report the results for the MAP ex- 
periment and compare them with the relevant predictions 
of Bouchet & Gispert (1999). A companion paper will de- 
tail the implementation and give results for the Planck 
Surveyor experiment that requires using a much more so- 
phisticated foreground model due to the higher resolution 
and higher sensitivity of that experiment. 

2. Full Sky Wiener Filtering 

Given a set of full sky maps observed in different frequency 
bands (at 22, 30, 40, 60 and 90 GHz for the MAP exper- 
iment), we want to recover the underlying CMB signal 
by combining the channels in an optimal way, that takes 
into account both the foreground contaminants and the 
instrumental noise. 

We write the signal as the sum of the true data and the 
noise in each band, d„ = s„ + n^, where each vector com- 
ponent is a pixel map. The filtered signal, s„, is a linear 
combination of the data, = W^d^, where W is the lin- 
ear filter used in the multi-frequency reconstruction. The 
Wiener filter is specifically obtained by minimizing the 
trace of the covariance matrix of the error map defined as 
We get the usual formula W = S{S + N)^ 1 
where S = (ss T ^ is the covariance matrix of the true sig- 
nal and N is the covariance matrix of the noise. The noise 
characteristics should be defined as precisely as possible 
for a given experiment, but the signal is a priori unknown. 
The prior probability given by the covariance matrix S is a 
key ingredient in the Wiener filtering method. We assume 
here that the signal in a given frequency band is a linear 
combination of several astrophysical processes including 
the CMB, the different Galactic emission processes, and 
other extragalactic foregrounds. This can be formulated 
as s„ = A vv ~x. v where x p is a reference template for a 
given astrophysical process. This assumption allows us to 
factorize the spectral and spatial properties of the signal. 
The observation matrix A can also account for the beam 
as a convolution operation on the input maps. The covari- 
ance matrix of the signal in the different frequency bands 
is readily obtained by S = ACA T where C is the covari- 
ance matrix of the templates (for uncorrelated processes, 
it is block diagonal). The last step in the component sep- 
aration is to recover from the filtered data the estimates 
of the astrophysical processes by a Least-Square fit to the 
recovered signal, namely 

-l , 



x p - (A T A)- 1 A T s v = CA T (AC A 1 



N) 



which is the formula for multi-frequency Wiener filtering 
( [Bouchet et al. 1996j ). 

One problem in analyzing full sky maps is the huge 
number of pixels that one must deal with. A direct pixel- 
based approach of the Wiener filtering technique is far 
beyond the capabilities of current and near-future super- 
computers. One possibility is to work with coefficients of a 
harmonic decomposition where convolutions (of the tem- 
plate by the optical beam) translate into mere multiplica- 
tions. Indeed previous simulation analyses all dealt with 



Fourier coefficients. For the full sky analysis presented 
here, we use instead spherical harmonics coefficients. This 
harmonic basis change is particularly useful if each as- 
trophysical process can be well approximated by a homo- 
geneous and isotropic random field with a given power 
spectrum. In this case, the Wiener filter can be computed 
for each multipole £ independently. 

Here we consider that, to first order, the Galactic emis- 
sion can be modelled as the sum of spatial templates for 
each astrophysical emission mechanism (i.e. dust, syn- 
chrotron) multiplied by their associated frequency depen- 
dence. At high galactic latitude this appears to b e a good 
approximation (see e.g. Bouchet fc Gispert 1999 , and ref- 
erences therein). At low galactic latitude, the situation is 
more complex as strong variations of the effective spectral 
index are observed in the Galactic disk. In addition, the 
presence of very bright sources at low latitude strongly vi- 
olates the assumption of Gaussianity that is required for 
Wiener filtering to be linearly optimal, while this assump- 
tion is more reasonable at high latitude. Moreover, the 
presence of these bright, highly localized features is likely 
to ruin any component separation analysis that assumes 
statistical isotropy. All these points, together with the fact 
that CMB emission is dominant only outside the Galactic 
plane, indicate that a careful removal of the Galactic plane 
region is necessary in the analysis of full sky CMB data. 

3. Working with a Galaxy Cut 

3.1. Building a New Orthonormal Basis on the Cut Sky 

The problem that arises as soon as one removes the 
Galactic plane region from a full sky map is the loss of 



orthonormality of the spherical harmonics basis (Gorsk 



1994). As was pointed out in the last section, a Wiener 



filter results from minimizing the trace of the covariance 
matrix E of the error map. If we write for a given map x 
that x = 7x, where Y are the spherical harmonics (Y; m 's) 
and x are the map's multipoles, we get 

Tr{E) = (e f e) = (e t Y T Ye) ^ (e f e) 

where e is the error map of any given process. Therefore, 
Wiener filtering the data in harmonic space is not equiv- 
alent anymore to Wiener filtering the data in real space. 
The usual way to solve this problem is to build a new or- 
thonormal basis on the cut sky. Gorski (1994) proposed a 
method for the COBE-DMR data based on a Cholesky de- 
composition of the coupling matrix Y T Y of the Y; m 's. In 
the course of our work, we found that this method works 
well for ^ max up to 40 - 50, but the Cholesky decomposi- 
tion failed to converge for higher multipoles because the 
coupling matrix becomes ill-conditioned due to numeri- 
cal truncation errors. We have therefore used instead an 
orthonormalisation scheme based on the Singular Value 
Decomposition of the Y matrix. The SVD can be written 
Y = UDV T where U T U = I, V T = V' 1 , and D is a 
diagonal matrix containing the singular values. This de- 
composition is numerically stable. The new basis is given 
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by the matrix U and the decomposition becomes x = C/ T x 
(if x can be represented on that basis, i.e. x = C/x, then 
(x T x) = (x T x) as desired). Vectors of the new basis that 
correspond to vanishing singular values can be dropped as 
their support is confined to the Galactic cut. Note that for 
high resolution maps, this new basis can be obtained at a 
reasonable computational cost (but only if the cut is sym- 
metric). We therefore restrict ourselves to a symmetric, 
isolatitude Galaxy cut (with \b\ > 20° for example). 



Input CMB map convolved with effective window 




w 



3.2. Power Spectrum Estimator on the Cut Sky 



One important ingredient of the Wiener filtering technique 
is the covariance matrix of the astrophysical processes. 
One possibility is to give an analytical form of the C(£) as 
a prior. In practice, however, it is more useful to compute 
the power spectrum directly on a first estimate of each as- 
trophysical t emplate (obtained for ins tance with an SVD 
method, see Bouchet fc Gispcrt 1999j ) and ultimately, do 
an iterative Wiener filtering by refining the estimate of 
the power spectra at each step. In this paper, we com- 
pute directly the power spectrum of the templates used in 
the simulated data, but by using only the pixels at high 
galactic latitudes. Tegmark (1995 ) described a method to 
obtain an optimal power spectrum estimator for each £ on 
an incomplete sphere by minimising the leakage between 



CMB error map 



multipoles (see Tegmark 1995, for details). We use this 
method to design "optimal masks" or "apodising func- 
tions" to compute the C(£) on our Galaxy cut templates. 
Our implementation was designed to work in the frame- 
work of the orthonormal basis on the cut sky described 
in the last section. An intermediate possibility would be 
to exploit the smoothness of the different power spectra. 
The co mputation of the p ower spectra in terms of band- 
powers ( Bond et al. 2000 ) might offer in this respect a 



good compromise between spectral resolution and compu- 
tational speed, but we did not test it. 



4. Numerical Techniques 

In this work, we use extensively the HEALPix pixelisation 
scheme^. We worked at a resolution n S id c = 256 with a 
pixel size of 13.7 arcmin^. This pixelisation choice is a 
little too coarse to extract all the CMB information from 
MAP since its 90 GHz channel has a 12.8 arcmin FWHM 
beam. As we shall see, however, it is enough to capture 
the essential features of the Wiener filtering possibilities 
for the MAP experiment. 

Due to the constraint that our Galaxy cut is symmet- 
ric with respect to the Galactic plane, we can address each 
m for odd and even £ separately. This rotational symme- 
try reduces the computational complexity from 0(Np ix ) 

down to 0(Np( x ). Computing the new basis represents 20 



http://www.eso. org~kgorski/healpix/ 




Fig. 1. Upper panel: The CMB input map, convolved with 
the effective window function of the experiment. Lower 
panel: The CMB error map. Note the different scale used 
here. The error is mostly white noise, except near the 
galactic cut. Several bright spots are also visible. 



minutes on a single processor of a SGI O2000 parallel sys- 
tem for n s ido = 256 and £ ma x = 512. Recomputing the 
new basis each time we need to analyze a map is not very 
efficient. We therefore store permanently the new coordi- 
nates, which requires approximately 1 GB space in double 
precision. Computing the "optimal masks" for power spec- 
trum estimates is much more computationally expensive. 
The core of the method is to solve a G eneral Eigenvalue 
Problem (Ax = XBx) (Tegmark 1995| ), and to compute 
only the first eigenvector. This process takes 45 minutes 
on the same parallel computer using 10 processors run- 
ning in parallel under MPI. The storage of these "opti- 
mal masks" also requires approximately 1 GB in double 
precision. Note that storage requirements are important 
here, since extrapolating our results to n s id c = 1024 and 
f max = 2048 leads to a total computational time of 48 
hours with 10 processors and a total memory requirement 
of 128 GB. We have, however, identified possible ways to 
decrease the CPU and memory requirements each by a 
factor of 4. 



Healpix divides the sphere in 12 main regions, each being 
hierarchically divided in n^ idc pixels 
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Fig. 2. Distribution of the CMB map reconstruction er- 
ror (dotted line), obtained by subtracting the recovered 
CMB map from the input CMB map smoothed with the 
effective window function for the CMB ("quality factor", 
see Bouchet & Gispert 1999) produced by Wiener filter- 
ing (see text). Also shown is the best-fit Gaussian (solid 
line) with its mean and rms. Note that the non-zero mean 
value shows the small level of contamination of the re- 
constructed CMB by the foregrounds monopoles. It will 



however not be an issue for the real observations since 



MAP is a differential experiment. 



Fig. 3. Comparison of the input sCDM angular power 
spectrum (upper thick solid line) to that of the recon- 
structed (Wiener-filtered, hence smoothed) CMB map 
(dots, middle noisy curve). The latter can be used to ob- 
tain an unbiased estimate of the input spectrum (thin 
solid line, on top of the input spectrum) by using the qual- 
ity factor produced by the analysis (see text). The lower 
(thick solid) curve shows the semi-analytical prediction of 
the spectrum of the reconstruction error by Bouchet & 
Gispert (1999| ), together with the actual estimate (thin 



solid) from our full-sky error map. 



5. Results and Discussion 



In this work we modelled the Galactic emission to be com- 
posed only of synchrotron and dust emission. Other forms 
of emission such as uncorrelated free-free emission or point 
and SZ sources were left out for simplification of the in- 
terpretation of the results. We will address those addi- 
tional components in a future paper. The dust emission 
was simulated using the Finkbeiner et al. DIRBE/IRAS 
composite 100 /im map (Finkbeiner et al. 1999), and the 
synchrotron emission using the Haslam 408 MHz map 
(Haslam et al. 1982) with spectral index —0.9. The dust 
was assumed to have a single 18 K component with emis- 
sivity index of 2, together with a correlated free- free like 
component of spectral index —0.15. 

Figure [l] shows the input CMB map, smoothed at the 
expected resolution of the reconstructed map, and the er- 
ror map, which is the difference of the input CMB map 
and the recovered one. As the latter shows, the overall 
reconstruction is excellent, although a few glitches associ- 
ated with low galactic latitude HII clouds are visible close 
to the Galaxy cut. This justifies a posteriori the removal 
of the Galactic plane from our analysis. Figure || gives the 
distribution of the reconstructed errors, whose bulk is well 
approximated by a Gaussian distribution with an rms of 
24 /iK. Note though the low level wings that are likely 
associated with the residual imprint of HII clouds. The 
corresponding pixels would probably have to be dropped 
out from any ensuing analysis. 



An important feature of Wiener filtering is that the 
quality factor(s) 3 of the instrument is an output of the 
analysis. The expected resolution of the reconstructed 
map is thus known directly and can be accounted for in 
later stages of the statistical analysis. Indeed we used it 
here to convolve our input CMB map to the expected res- 
olution of the recovered map. 

For cosmological purposes, the most important statis- 
tic is the angular power spectrum of the CMB. We there- 
fore show the power spectrum of the recovered CMB maps 
together with the spectrum of the reconstruction error in 
figure [|. The spectra are computed for each i using the 
minimal leakage apodising masks on the incomplete sphere 
(see section 3.2). The curves abruptly end at i = 512, the 
resolution limit of the present analysis, where incidentally 
the error is of the order of the CMB signal. Note however 
that a box car averaging of the C(£) with A£ around 10 
will reduce the noise level, and therefore a strong (and thus 
useful) cosmological signal is still expected above I — 512, 
likely up to £ — 1024. We will extend this work in a com- 
panion paper to a spatial resolution of at least i = 1024. 
The main conclusion of the current work is that we con- 



firm the semi-analytical predictions of Bouchet & Gispert 



| 3 [ The quality factors are the generalized instrumental win- 
dow functions for each recovered astrophysical process, taking 
into account the angular resolution and the detector noise of 
eve ry channel, as well as the contamination of other processes, 
see Bouchet & Gispert 199£ . 
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(1999), but using this time a full sky analysis of simulated 
MAP data. 
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